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An investigation was conducted to determine the 

feasibility of obtaining the short period and phugcid 

stability derivatives from one maneuver, 

"simultaneously”. It was concluded that the maximum 

likelihood identification program SCIDNT I showed 

great promise in obtaining the short period and 

phugoid stability derivatives from one 

maneuver. Extraction of the short period stability 

parameters in the presence of the phugoid was 

easy, straight forward, and yielded results similar to 

those obtained from pure short period data. Estimation 

of the phugoid stability parameters was possible when 

they were estimated in conjunction with 

Z , X , a , and 0 . It was recommended that a new set 
o o o o 

of data be obtained at a flight condition where the 
phugoid is at least moderately damped and that this 
data be analyzed to resolve the present anomalies. 
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I. INTRODUCTION 



A. BACKGROUND 

Flight testing of military aircraft has always been a 
time consuming and expensive proposition. The relatively 
recent introduction of the new Flying Qualities 
Specification (MIL-F-8785B) has introduced many new 
parametric requirements which require that many • of the 
aircraft airframe parameters be determined to insure 
specification compliance. As a result of these new 
reguirements the Naval Air Test Center (NATC) initiated a 
research program to develop airframe parameter 
identification technology for use in flight testing Navy 
aircraft. This program was initiated in 1971 between the 
Naval Air Systems Command (NAVAIR) and NATC. At the same 
time the Office of Naval Research (ONR) issued a contract to 
Systems Control, Inc. (SCI) to advance the state-of-the-art 
in parameter identification. In 1973 a concentrated effort 
was begun by the Navy to continue the development of this 
new technology by forming a joint program between 
NATC, NAVAIR, ONR and SCI. As a result of this development a 
digital FORTRAN computer program, SCIDNT-I (hereafter 
refered to as SCIDNT) was written for NATC to extract 
aircraft stability and control coefficients from flight test 
data. Ref. 1. It was hoped that by using this computer 
program it would be possible to significantly reduce the 
flight test time in the air, and the ground data reduction. 
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B. PURPOSE 



The purpose of this research was to develop parameter 
identification test methods for determining aircraft phugoid 
and short period characteristics from one aircraft 
maneuver. This research further entailed the determination 
of proper aircraft maneuvers/pilot inputs and detailed data 
analysis proceedures using the advanced maximum likelihood 
identification technique. 



C. DESCRIPTION OF AIRPLANE 

The test airplane was demonstration F-I^JA ship number 8 
BU. NO. 157987. The airplane was essentially typical of a 
production F-14A airplane relative to these tests. An 
aircraft description is contained in Appendix A. Photographs 
and a three view drawing of the test airplane are shown in 
Appendix A , figures 1, 2 and 3. 

D. KAXIIiUM LIKELIHOOD PARAMETER IDENTIFICATION THEORY 

A widely used rule for choosing the sample statistic to 
be used in parameter estimation is called the method of 
maximum likelihood Ref. 2. This method, in addition to 
prescribing the statistic wh.ich should be used, provides an 
approximation of the distribution of that statistic so that 
approximate confidence intervals may be constructed. The 
method of maximum likelihood makes use of the 
sample- likelihood function. Assume we observe or record a 
set of observations on a system, whose model has p unknown 
parameters 6. -We then can write the joint probability 
distribution of a random sample for any given set of values 
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of the parameters 0 from a feasible set 0. We then can 
assign a probability p(z|0) to each outcome z. If we 
observe or measure the outcome of an actual experiment to be 
z we then would like to know which set or sets of values of 
0 might have resulted in these observations. This concept is 
contained in the likelihood function ac (0|z). This function 
is of fundamental importance in estimation theory because of 
the likelihood principle of Fisher and others which states 
that if the system model is correct, all information about 
unknown parameters is contained in the likelihood 
function. The maximum likelihood estimator of 0 is the value 
0 which causes the likelihood function oC (0|z) to be a 
maximum [2]. 



% = max (6|z) 
0£d 



From a purely mechanical viewpoint this method can be 
summarized as follows: 



Find the probability density functions of the 
observations or measured values for all possible 
combinations of the unknown parameter 
values. Select the density function whose value is 
highest among all density functions at the 
measured values of the observations. The 
corresponding parameter values are the maximum 
likelihood estimates. 



Assume 0 can take three possible values: 9,0, and 

12 

0^. Further assume that the probability density functions of 

observations z for these three values of 0 are as shown in 

Figure 1. Then, if the actual observation is z, 0 is the 

2 

maximum likelihood estimate of 0. 
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z 

observed 

Figure 1. fiaximum Likelihood Estimates 

Owing to the aonotonic, one-to-one relationship between 
the likelihood function and its logarithm, i.e. they have a 
maximum at the same value 0, it is more convient to work 
with the negative of the logarithm of the likelihood 
function. Starting from a priori values, the parameters are 
updated so that the value of the negative log- likelihood 
function J at the observed value of the outputs decreases 
continuously . 

As sample statistics or functions of the 
observations, the estimators are, before the experiment, 
random variables and can be studied as such. In fact, the 
major advantage of maximum likelihood estimators is that 
their properties have been thoroughly studied and are well 
known, virtually independently of the particular density 
function under consideration. Additionally, the maximum 
likelihood method can be used with linear or nonlinear 
models in the presence of process and measurement noise. The 
estimator of 0 is approximately normally distributed for 
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unbiased, efficient. 



large n, and is asyia ptotically 
consistent, and sufficient, [2]. 

The details of the maximum likelihood identification 
proceedure are presented in the following section. 



E. MAXIMUM IIKELIHOOD IDENTIFICATION OF PARAMETERS 

The derivation follows closely that presented in Ref. 
3, and is presented herein for completeness and reader 
convenience. The initial equations will be concerned with 
general ncnlinear systems. 

Consider the general nonlinear aircraft equations of 
motion 



X = f(x, u, e, t) + P(6, t) w 0 < t < T ( 1 ) 



E(X(0)) = X (0) 

o 



E{(X(0) - X (0)) (X(0) - X (0)) } = P (9) 



( 2 ) 



where 



X (t) 


is 


n 


X 


1 


state vector 


u(t) 


is 


1 


X 


1 


input vector 


e 


is 


P 


X 


1 


vector of unknown parameters 




is 


n 


X 


g 


process noise distribution matrix 


w(t) 


is 




X 


1 


random process noise vector 



Sets of m measurements y(t ) are taken at discrete times t 

k k 
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y(t ) = h(x(t ), u(t ), e, t ) + v(t ) ( 3 ) 

h k k k k 

k = 1, 2, 3,. ... ,N 



w’(t) and V (t ) are Gaussian random noises [2], with the 



following properties 



B(w(t)) = 0 E(v(t )) = 0 

k 

T 

E(w(t) V (t )) =0 
k 



T 

E(w(t) w (^')) = Q (0,t) (t-T ) 



T 

S(v(t .) V {t )) = R(0,t .) <? ., ( 4 ) 

3 k 3 3k 

The unknown parameters are assumed to occur in the 

functions f and h and in matrices P, Q, R, P ^ and x . In 

' o o 

the following analysis, the model and the functional form of 



f and h is assumed known correctly. 

The set of observations y (t ),y{t ),...., y(t ) 

1 2 N 

constitutes the outcome Z in this case. The likelihood 



function for 0, which has the same form as the probability 
of the outcome z for a certain value of parameters 0, is 
given by 



(6iz) = p(zi0) 

= P(y<tj , y(t^) ,...,y{t^) |0) 
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( 5 ) 



where 



^ (0|z) = p(Y |0) 

N 



Y — {y('fc )r*«*/y(^ )}# 2, 

k I k 



p(Y^|9) = p(y(y P(V^.,ie) 






which by successive application of Bayes rule, [2], becomes 



N 

P(Y |0) = TT p(y (t.) |Y , 0 ) ( 6 ) 

N ' ' 11-1 

i=1 



The log-likelihood function is 



N 



log (^(01z)) = ^ log 

i-1 



{P(y(t ) |Y / 6) } + constant 
i i-1 



( 7 ) 



To rind the 
0, the mean 



probability distribution of y{t ) given Y 

i i-1 

value and covariance are determined first. 



and 



E(y(t.) I Y. , 0) 6 y(i|i-l) ( 8 ) 

1 1-1 

The expected value or the mean is the best possible estimate 
of measurements at a point given the measurements up to and 
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including the previous point. 



cov(y(t.)|y. ,0) = E{(y(t.) - y (i|i-l)) (y(t.) - y (i| i-1) ) ’’} 

11-1 1 1 



AE{v'(i) V^(i)’’} 

A B(i) ( 9 ) 



(i) are called the innovations at point i and B(i) is the 
innovations covariance. Since 



y (t.) - E(y(t.) |Y. ,0) = \^(i) ( 10 ) 

1 11-1 

it follows that y (t } given Y and 0 have the same 

i i- 1 

distribution as v' (i) . It has been shown that as the 

sampling rate is increased, the innovations V (i) tend 
towards having a Gaussian density. Assuming a sufficiently 
high sampling rate, the distribution of V" (i) and, 
therefore, y(t_) given Y^ and 0 is Gaussian, i.e.. 



p(y(t.)iY, ,0) 
1 1-1 



r 1 ^ T 

exp - — ( 1 ) B-i ( 1 ) 

ra /2 1/2 

(2 TO |B(i)| 




( 11 ) 
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log {p(y (t.) !Y. ,0)} 

1 x-1 




(1) B-Ml) V(l) 



1 

-- loglB(i) I 
2 



+ constant 



( 12 ) 



The log- likelihood function of equation ( 7 ) can be written 
as 



N 

log ( ^( 0 | 2 )) = - y ^ {/ (i) B-i(i) v' (i) 

i=1 

+ logIB (i) 1} { 13 ) 

An estimate of the unknown parameters is obtained by 
maximizing the likelihood function or the log- likelihood 
function from the feasible set of parameter values. 



0 = max log ( ^ (0 | z) ) ( 1^ ) 

©£0 

r 1 N 

= max - y ^ iV’ (i) B~Mi) (i) + log I B (i) 1} 
0£-0 ^ i=1 •' 



( 15 ) 



The log- likelihood function depends on the innovations 
and their covariance. To optimize the likelihood function, a 
way must be found for determining these quantities. Both 
innovations and their covariance are outputs of an extended 
Kalman filter. 
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The extended Kalman filter is conventionally divided 
into two parts. In the first part, called the prediction 
equations, the state equations and state estimate covariance 
equations are propagated in time from one measurement point 
to the next. In the second part, called the measurement 
update equations, the measurements and associated 
measurement noise covariances are used to improve state 
estimates. The covariance matrix is also updated at this 
point to reflect the additional information obtained from 
the measurements. 

1 • P re d iction Eq uat i on s 

The state prediction is done using the equations of 

motion. Appendix B. Starting at time t with current 

i^1 

estimate x(i-1|i-1) of the state x (t ) and the covariance 

i 

P(i-lli-l), the following equations are used to find the 

predicted state x(i|i-1) and the associated covariance 
P(i|i-1) . 



d 

— x(t|t. ) = f(x(t|t. ),u(t),0,t) ( 16 ) 

dt 1-1 1-1 



P(t|t. ) = F(t) P{tlt. ) + P(t|t. ) F^(t) 

1-1 1-1 1-1 






t < t < t 
i-1 i 
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The n X n matrix F is obtained by linearizing f 
about the best current estimate 



F(t) 



3f(x(t|t. ),u(t),e,t) 

X- 1 

ax(tit. ) 

x-' 1 



( 18 ) 



Using ( 16 )-( 18 )/ we can obtain 



X (t It . ) X (il i-1) 

X X- 1 



and 



P (t. It. ) A P(i| i-1) 

X X- 1 



( 19 ) 



Thereafter, the measurement update equations are used. 

2- Measurement Update E quations 

The covariance and state estimate are updated using 
the measurements. The necessary relations are derived in 
Ref. 4, and .are presented here without proof. The innovation 
and its covariance are 



V"(i) = y(i) - h(x (i|i-1) ,u(t.) ,0,t.) ( 20 ) 

X X 



T 

B(i) = H (i) P(i|i-1) H (i) + R ( 21 ) 
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where H is obtained by linearizing h. 



H(i) 



3h(x(i|i-1) ,u(t. ,e,t) 

1 1 

3x (il i-1) 



( 22 ) 



The Kalman gain and the state update equations are 



T 

K(i) = P (ili-1) H (i) B-i (i) 



( 23 ) 



x(i1i) = x(i|i-1) + K(i) v" (i) 



( 24 ) 



Finally, F(iji), the covariance of error in updated 
state, is obtained by 



P(i|i) = (I-K(i) H(i)) P(i|i-1) 



( 25 ) 



One can now return to the time update or prediction 
equations. 

Optim iza tion Procedure 

Many possible numerical procedures can be used for 
this optimization problem. Modified Newton-Raphson or 
Quasilinearization have been found by experience to give 
quicker convergence than most procedures such as the 
conjugate gradient or the Davison method. The modified 
Newton-Raphson method is a second order gradient procedure 
requiring computation of first and second order partials of 
the log- likelihood function. 
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j ,k = 1,2, ... ,p 



The last three terms in the equation for the second 
partial of the log- likelihood function involve second 
partials of innovation and its covariance. These terms are 
usually dropped which yields the second partial 
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approximation 
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From equations 20 and 21 we can obtain the gradients 

of innovation and its covarience for parameter 0 as, 

3 
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Recuisive equations can be obtained for gradients of 
the predicted state and its covariance. This is done in 
stages by using the prediction and measurement update 
equations of the Kalman filter. Differentiating equations 
16-18 with respect to 0 yields 
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j = 1,2,. - . ,p 



The sensitivity functions, used in the SCIDNT integration 

routines, are updated at measurement points by 

differentiating equations 22 - 25 with respect to 0 . 
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The negative of the matrix of second partials of tne 
log-likelihood function is called the information matrix 
M. The step size A 0 for parameter estimates is given by 
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The information matrix provides a lower bound on parameter 
estimate covariances, i.e., 

T 

E (e-e) ( 0 -e) > ( 39 ) 



This is the Cranter-fiao lower bound. The maximum likelihood 
estimates approach this bound asymptotically. 

h • iinear Syst ems 

In a linear system, the functions f and h are 
defined as 



f(x,u,0,t) ^F(0,t)x *■ G(0,t)u 



and h(x,u,0,t) ^H(0,t)x + D(0,t)u ( 40 ) 

For the linear case, it has been shown that if the 
model is correct, the innovations are white and have 
Gaussian density at the true values of the parameters. The 
assumption of fast sampling rate is, therefore, not 
necessary. 

The basic algorithm for solution of the linear case, 

is the same, however, some of the equations can now be 

simplified. The equivalence between ?(t) and H (i) of ( 18 ) 

and ( 22 ) and F(0,t) and H{0,t ) of ( 40 ) is 

i 

obvious. Equations ( 16 ) and ( 20 ) now become 
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and 
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Equations ( 29 ) and ( 31 ) can be written as 
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All Other equations remain the same. 

There is a considerable reduction in the computation 
requirement for time-invariant linear systems. In this 
case, matrices F, G, H, D, p , Q and R and their derivatives 
with respect to the parameters are constant. 
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5 • Ti rae Invaria nt Linear Systems In Statis tical Stea dy 
State 

In many aircraft applications, the Kalman filter is 
in steady state for the duration of the experiment. This 
occurs when the Kalman filter is in operation for a 
sufficiently long time and the process and measurement noise 
covariances are constant. The Kalman gain and the 
innovations and the state covariances approach constant 
values. The time update and measurement update equations for 
the covariances are 



4-p(tit. ) = pp(tit. ) + p(tit. 

dt 1-1 1-1 1-1 ' ' 

( 45 ) 
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By definition of the steady state 



P (i-1 1 i-1) = P(i|i) 



Therefore, from ( 45 ) 
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Using ( 46 ) , ( 47 ) and ( 49 ) , we can solve for P{i|i-1) 
and then find K and B. Also, it can be shown that 
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Thus, it is possible to solve for 
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is a linear equation in 
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and the coefficient of the 



unknown matrix does not depend on the parameter 
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0 . Thus, the sensitivity of the state covariance matrix can 

j 

be determined very quickly for all parameters. Once the 
sensitivity of P, K and B for unknown parameters is 
determined, only state sensitivity equations need to be 
updated. The computation of state sensitivity functions can 
be reduced to many fewer equations. 

An approximation simplifies the problem further. The 
unknown parameters are defined to include elements in K and 
B matrices instead of Q and R. Optimizing the log-likelihood 
function for parameters in B gives 



N 

B = — ) V'(i) / (i> 

i=l 



( 54 ) 



The gradient of the log- likelihood function with respect to 
other unknown parameters is 
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The sensitivity of innovations to parameters is determined 
using the following recursive equations 
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( 57 ) 



a v"(i) 
as. 



3H , 3x(i|i-1) 3D 

x{i|i-1) - 



ae . 

3 



ae. 

3 



ae. 

3 



( 58 ) 



3x (i|i) 

~w. 

3 



3x(i|i-1) 3K a V"(i) 

+ \/(i) + K 



ae. 

3 



ae. 

3 



ae. 

3 



3 — 1#'2/..»/p 



Note that 



=0 if 9 is not an element of K matrix 

ae . 3 

3 

( 59 ) 

= I if 9 A K 

j'k* j j'k* 



39 



where I is a matrix of all zeroes except a 1 at the 

j'k* 

j', k' position. 

This approximation simplifies the optimization 
considerably. However, this usually leads to an 

overparameterized model. In other words, once K and B are 
determined, it is not possible to find any corresponding 
p, Q and S which have the desired (known a priori) 

structure. A good estimate of elements in P r Q and R 
matrices can be obtained by a least-squares type 

approach. The fit to the observed data is better than with 
true values of parameters but the parameter estimates do not 
have minimum variance. This approximation is not good in 
aircraft application where the structure of Q and R is known 
fairly well, but is excellent where there are many process 

noise sources and the characteristics of both Q and R are 

relatively unknown (e.g., economic systems). 

6* Lik elihood Wirh No Proce ss Or Measu r eme nt 

Noise 

The maximum likelihood method can be simplified when 
either process noise or measurement noise are 
absent. Assuming no process noise we can proceed as 
follows. 

If the process noise is zero and initial states are 
known perfectly, i.e. , w(t) and P(0) are zero, the 
covariance of the error in the predicted state is also 
zero. It is clear from ( 23 ) that the Kalman gains are 
zero. The innovations are the output error, i.e.. 



V"(i) = y(i) - h(x(tj ,u(t.) ,e,tj ( 60 ) 

111 
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and the innovation covariance is ( 21 ) 



B(i) = R 



( 61 ) 



the log- likelihood function is. 
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log ( 5 e(e|z)) = ^ R-^V(i) + log |R| 
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( 62 ) 



which on optimizing for unknown parameters in R yields 



N 

R = — y \/(i) V' (i) ( 63 ) 

N ^ 

i=1 



The equality in ( 63 ) holds only for those elements of R 
which are not known a priori. For instance, even if R is 
known to be diagonal, the right hand side matrix will not be 
diagonal in general, but the off-diagonal terms should be 
ignored before they are equated to R. Using ( 63 ) in ( 62 ) 
we obtain 



N 

log (oC(0|z)) = -yL ^ 

i=1 

+ constant ( 64 ) 
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The optimizing function is the same as that for the output 
error method except that the measurement noise covariance 
matrix is determined using ( 63 ) and is used as the 
weighting matrix in the criterion function. In the output 
error method, the measurement noise is assumed known and the 
weighting function is arbitrary. 

The first and second derivatives of the 
log-likelihood function with respect to unknown parameters 
are 
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The terms in the second derivative are approximated as 
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With no measurement noise we can proceed as 
follows. If all states are measured with no noise, the 
covariance of the error in state estimates is zero at the 
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beginning of any time update. 



P (i-1|i-1) = 0 



and 



x(i-1|i-1) = x(i-1) 



( 68 ) 



It is easy to show in this case that for fast sampling the 
log-likelihood function is quadratic in the difference 
between measure values of x and f(x,u,6,t). The method 
reduces to the equation error method, the weight W being 
chosen as 



w 




f(x,u,e,t)) (X 



T 

f(x,u,e,t)) dt 



( 69 ) 



Thus, the maximum likelihood method and equation-error 
methods are equivalent except for the technique for choosing 
the weighting matrix. 



F. AIRCSAFT PARAMETER IDENTIFICATION USING THE MAXIMUM 
LIKELIHOOD METHOD 

In applying the maximum likelihood theory it is 

necessary tc consider the aircraft equations of motion in 

discrete form. It was assumed that the structure of the 

model was known. The vector of unknown parameters from x (0) , 
0 

0 , G, P , H , Q, and R was denoted by 6. The 



as 



problem, therefore, was to identify 0 using the maximum 
likelihood method. In constructing the problem it was 
further assumed that the aircraft which had been flying for 
seme time was in a steady state condition, and process noise 
(wind gust effects) was not present. 

In the computer implementation, the negative 
log-likelihood function was minimized using a Gauss-Newton 
gradient procedure. A least squares routine was utilized to 
provide approximate airframe parameter start-up 
values. Additionally the results of the SCIDNT computer 
program were utilized to obtain frequency and damping. 



G. SCOPE OF RESEARCH 

Research was limited to the longitudinal case since both 
available computer and research time were limited. Although 
additional data were desired, the project airplane was 
substantially damaged by fire and was removed from further 
test activity. Further data gathering with the F-14A 
airplane was, therefore, not possible. Specific 
longitudinal stability derivatives identified are given in 
the Results and Discussion section. There were no airplane 
flight restrictions which affected the test methods or 
limited the scope of research. The research was limited to 
establishing base line data for the short period and phugoid 
modes and thence to determine the feasibility of identifying 
these modes simultaneously. 



H. HETHOD OF TESTS 

In general test proceedures were in accordance with 
those contained in the U. S. Naval Test Pilot School Fixed 
Wing Stability and Control Manual, Ref 5. There 
were, however, some non-standard test techniques utilized in 
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exciting the desired airplane modes. The short period mode 
was excited by the standard doublet method and additionally 
by applying a sinusoidal control input at varying 
frequencies. The phugoid mode was excited by the Delta 
Airspeed method in all cases, [5]. 

In exciting the phugoid and short period modes 
simultaneoulsy, a combination of standard test techniques 
were utilized. The "simultaneous" technique consisted of 
establishing the airplane in a phugoid maneuver and daring 
the first quarter cycle of the maneuver applying a doublet 
input to excite the short period mode. 

The data system utilized within the airplane was 
specifically designed to be compatible with the Naval Air 
Test Center Real-Time Telemetry Processing System (RTFS) . A 
brief description of the airplane and ground data systems is 
presented in Appendix C. Data reduction was performed on the 
RTFS Xerox Sigma 9 computer in the batch processing 
mode. Raw airplane data were converted to an Engineering 
Unit tape and used as input to the digital parameter 
identification program (SCIDNT) . In addition a least squares 
parameter identification routine was utilized to generate 
parameter start up values for SCIDNT. This least squares 
routine extracted a set of approximate stability derivatives 
by minimizing the error between the model and actual 
airplane state vector time histories. 
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II 



• results and discussi on 



A. SHORT PERIOD AIRFRAHE PARAMETER IDENTIFICATION 
Ai r plane Test Cond iti ons 

Initial trimmed test conditions for each control 
input are presented in Table I. Prior to each control input 
the airplane was trimmed for steady state "hands off" flight 
[5]. 



Altitude 
(ft HPC) 

Sine Have Input 29 925 

Doublet Input 30 503 



MACH CG Loading 

(% MAC) 

( 1 ) 

.6149 11.18 B 

.6289 10.985 B 



(1) . Four Sparrow III missiles (AIM-7) on stations 3, 4, 5 

and 6 



Table I. Short Period Test Conditions 



2. Metho d Of Investigat ion 

The airplane model equations used were those shown 

in Appendix B. In order to isolate the short period mode 

from the longitudinal equations of motion the state vector 

parameters oc. , q, and 0 were made insensitive to u. Thi-s was 

accomplished simply by fixing 3^ , X , X , Z , M and X^ 

u q u u e 
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equal to zero 


and not identifying any 


of them. 


This 


equation 


modification 


resulted in 


the 


clasical 


short period 


approximation 


Ref. 6. 










To provide SCIDNT 


with 


approximate 


airframe 


parameter start-up values, the 


least 


squares 


routine was 


utilized. For 


the short period 


case 


the least 


squares start 


up values were 


very good, a comparison of a typical 


set of 


least squares 


parameters and 


the 


final SCIDNT 


parameter 


values is shown in Table II. 













Z 


Z 


H 


M 


2^ 






CC 


g 


CC. 


g 


<^e 




Least 


-.4909 


0. 9624 


1 

• 

00 


-.9085 


0.210 


6.09 


Squares 

SCIDNT 


-.4903 


0. 9805 


-1.95 


-.795 


0.278 


5. 542 



Table II. Comparison of Least Squares Parameter 
Values and SCIDNT Parameter Values 

In actual practice, however, airplane data was 
recorded in long bursts which usually included three or more 
maneuvers , executed with near identical initial 
conditions. As a result once a set of start up values were 
determined for a given test condition it was unnessary to 
determine a set of start up values for other like 
maneuvers. 

The actual parameter identification proceedure , for 
the short period mode, was relatively simple once 
familiarity with the digital program was attained. Upon 
obtaining a set of "reasonable" initial parameter 
estimates, all six short period parameters were estimated in 
an iterative fashion. Here "reasonable" denotes proper 
parameter sign and magnitude. Each computer run was made 
with the number of iterations and step cuts both initialized 
at ten. During a given iteration the program attempted to 
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minimize the likelihood function. If the likelihood function 
increased/ a series of step cuts would be taken until the 
likelihood function again decreased or the maximum number of 
step cuts were exceeded. The program would then continue to 
iterate until the desired number of iterations/step cuts 
were reached. Due to program logic one more step cut than 
requested was executed. For the short period 
mode, however, convergence was obtained early in program 
execution. In general, iterations conducted with doublet 
control input data converged very quickly and began 
executing step cuts until the allowable number was 
exceeded. Figure 2 shows the change in the likelihood 
function during iteration for a doublet input. Program 
convergence had been essentially reached by the fourth 
iteration at which time a series of step cuts were made. 



Iteration 




Doublet Input 



The convergence characteristics using sinusoidal 
input data were somewhat different. During the program 
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iteration the likelihood function continuously decreased and 
no step cuts were taken. Figure 3. Program convergence for 
this case was evidenced by a leveling out of the likelihood 
function, but not until the eighth iteration. It was felt 
that since the data length for the sinusoidal input was 
longer than the short period data length the larger number 
of iterations until program convergence was not abnormal. 



I teration 




Figure 3. Likelihood Function Change During Iteration For 

Sinusoidal Input 



After convergence was obtained for the parameter 
estimates, the Angle Of Attack Scale Factor (K^ ) was 
determined. was a scaling factor in the angle of attack 
measurement equation. Appendix B. In estimating two 

terms in the model were actually affected, itself and 

the variable (1 K /U) which was used to multiply q in the 

C\, 

measurement equation. Since, in essence, two terms were 

being estimated a modified form of parameter estimation was 
utilized. In the term (1 . K /U) the value of was 

initially set to 1.0 and the single term was 



determined. This new value of K , was then substituted into 

oc 

the term (1^ K /U) and again a set of iterations were 
performed to determine an updated . This method of 

iteration was continued until there was essentially no 
change in , at which time the value was said to have 

converged. In most cases only two iterations were reguired 
to obtain a converged K^. . 

Upon determining a value for ,a final set of 

iterations were • performed to determine values for 

Z,X,M,e.At this point, the short period mode time 
o o o o 

history fits were nearly identical and the estimate 

statistics indicated a high degree of confidence. Attempts 
at further iterations to improve the estimates resulted in 
using a lot more computer time , with very little 

improvement. Therefore, the procedure utilized to estimate 
the short period mode was to first estimate the short period 
parameters, second estimate , and finally, estimate the 

initial values. 

3 . Data An alv si s 

Data from typical short period estimation runs for a 
doublet and sinusoidal input are presented in Tables III and 
IV. Typical short period time history fits are presented in 
Appendix D Figures 1 and 2, for the doublet and sinusoidal 
inputs respectfully. It should be noted that the scales for 
the short period time histories are different for the 
doublet and sinusoidal inputs. Furthermore, since the 
estimated stability derivatives were in dimensional form the 
actual derivative values differed slightly from maneuver to 
maneuver . 

The data obtained wit-h the doublet input in most all 
cases converged to believable derivative estimates, i.e., 
expected sign and magnitude, with favorable iteration 
statistics. For example, refering to Table III, the value of 
was -1.844 with a standard deviation of 0.0061 and a 
final step size of -0.0146 which is significantly lower in 
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magnitude than the estimated parameter. From probability 

theory we can say that the probability is 99.74% that the 

value of is - 1.844 ±0.0183 which is well within 

cc 

acceptable engineering tolerances. This same analysis was 
carried out for the other estimated derivatives with similar 
results. 



Derivatives 
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oC 


Z 
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M 

cc 


M 

q 




"-e 


Estimate 


-.4138 


1.099 


- 1.844 


-.82 6 


0. 1628 


5.098 


Standard 

Deviation 


0.0024 


0. 0038 


0. 0051 


0.0037 


0.0055 


0.0104 


Parameter 
Step Size 


-.0006 


-.012 


-.0146 


-.0001 


0.0143 


0.028 


Eigenvalue 

Spread 






1.28 


X 10 * 






Likelihood 

Function 

Value 






- 5 . 1 


X 103 






Time 

History 

Matches 


cr = 0. 

ec 


00396 


cr = 

u 


0.564 


cr = 0. 
g 


0039 






(T = 
e 


0.0043 


67- = 2. 
n 

z 


33 





Table III. SCIDNT Estimates of Longitudinal Short Period 
Stability Derivatives Using a Doublet Input 



The eigenvalue spread, that is, the difference 
between the largest eigenvalue and the smallest eigenvalue 
of the information matrix, was 1.28 x 10* which was less 
than the maximum spread allowed by SCIDNT (lO'O). All 
parameters converged indicating no identif iabilit y 
problems, and the eigenvalue spread was therefore not 
considered significant. The likelihood function value was a 
large negative number which converged rapidly (Figure 
2) , further indicating no identifiability problems. 
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As a final check of stability derivative correctness 
a series of time history fits were constructed. In 
constructing these plots both the model and airplane were 
excited with the same control input and the airplane time 
history plots superimposed on the same scale. The time 
history fit for the doublet input was, in general, very 
good, with the difference between the measured value and 
estimated value within the instrumentation measurement error 
tolerances Ref. 7. 

Referring to Table IV, it can be seen that typically 
the estimated derivatives using a sinusoidal input converged 
to appropriate (sign and magnitude) stability derivative 
values. 
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M Z^ 
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Estimate 


-.4303 


1. 1229 


-2. 007 


-.7615 0.1532 


5.2274 


Standard 

Deviation 


0.0042 


0.0061 


0.013 


0.007 0.0045 


0.0131 


Parameter 
Step Size 


0.0007 


-.0016 


-.0043 


-.0061 -.0037 


0.0332 


Eigenvalue 

Spread 






1.67 


X 105 




Likelihood 

Function 

Value 






-7.323 x103 




Time 

History 

Matches 


• 
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cr = 
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• 

II 
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006 






0~ = 


0.01 


(F = 2.32 
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Table IV. SCIDNT Estimates of Longitudinal Short Period 
Stability Deravitives Using A Sinusoidal Input 



In general, however, the parameter standard deviations and 
step sizes were not as favorable as those obtained with the 
doublet input for these two examples. It was 
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possible, however, to select another pair of maneuvers and 
reach the opposite conclusion, ie. sinusoidal input more 
favorable. Since all parameter estimates for both the 
doublet and sinusoidal inputs converged it was concluded 
that the slight difference in statistical data was not 
significant. If a determination as to which control input to 
use were required, it was felt that the primary 
consideration would be ease of pilot input. If an autopilot 
device were available which could input the sinusoid then 
this input should be used, however, with no autopilot system 
installed the doublet input would provide excellent 
results . 

A consistent problem that occured throughout the 
research was concerned with the initial conditions prior to 
the control input. If the initial condition of the airplane 
was not stabilized prior to the control input or if a data 
burst was taken in the middle of a control input, no time 
history match was possible. According to theory it is 
possible to extract the airframe parameters at any time 
during a maneuver, this, however did not appear to be 
true. Figure 3 Appendix D shows the results of estimating 
the parameters starting from non-zero initial 
conditions. This time history fit was generated after all 
iterarions on both parameters and initial conditions were 
completed. The parameter estimates did not converge and the 
most powerful derivative, H , was much larger in magnitude 

CCr 

than any ether estimates with a value of -3.7. These results 
occured any time an attempt was made to start SCIDNT from an 
initial condition other than steady state. It was concluded 
that, within the scope of this research, meaningful 
parameter estimates were only obtainable when started from 
an initial steady state condition. 

A comparison of short period charac teristics 
obtained from Navy Preliminary Evaluation (NPE) data and 
those determined by SCIDNT is presented in Table V. 
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Description 


Damping 

Ratio 


Frequency 




^sp 


CD 

n 

sp 

(rad/sec) 


NPE 


0.35 


1.49 


SCIDNT-Doublet 


0.403 


1.54 


SCIDNT-Sinusoidal 


0.371 


1-61 



Table V. Comparison of Longitudinal 
Short Period Characterist ics 



It was concluded that, within the scope of this 
research, the parameter identification program SCIDNT was 
easy and straightforward to use in determining the short 
period mode. In addition, it was concluded that it was not 
possible to obtain accurate estimates when starting from 
non-steady state initial conditions. 



B. PHUGOID AIRFRAME PARAMETER IDENTIFICATION 

1 ♦ Ai r plane T est C ond i tion s 

Initial trimmed test conditions for the phugoid 
maneuver are presented in Table VI. Prior to the control 
input the airplane was trimmed for steady state ’’hands off” 
flight [5]. 



Phugoid 



Altitude MACH 
(Ft KPC) 

Input 30 055 0.614 

Table VI. Phugoid Test 



CG 

(% MAC) 
10.94 

Conditions 



Loading 

B 
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2* Me th od Of In vestig ati on 

The airplane model equations used were those shown 
in Appendix B. The derivitives and X^ were fixed at 

zero values since their effect on the phugoid was 

small. Additionally, it was felt, fixing these terms would 

facilitate the identification of the dominant phugoid 

terms. Finally from Appendix B it can be seen that the terra 

X^ was defined as (X - oc U ) and therefore due to its 
g .goo 

relatively large size was approximated by cC U and 

o o 

fixed. The parameters which were identified during the pure 

phugoid maneuver were therfore, Z , X , and M . 

u u u 

Since there was no short period input during the 
pure phugoid maneuver it was not possible to effectively 
identify the short period parameters. Additionally the 
program solution diverged without a set of predefined short 
period parameters and it was necessary to obtain parameters 
from a prior run. A set of short period parameters 
were, therefore, obtained from a previous identification run 
and fixed throughout the phugoid parameter identification 
runs. 

In an attempt to provide SCIDNT with approximate 

airframe parameter start-up values the least squares routine 

was utilized. For the phugoid case the least squares start 

up values were unusable. For example, least squares values 

for X was -2.6697 and was -2h0.5, both of which are 

cc <^6 

terras of small absolute value. As an additional check, the 

least squares parameter estimates were used to generate a 
set of time histories. For these time hisxorys the standard 
deviation of airspeed and normal acceleration was 49.12 
ft/sec and 5.44 ft/sec^ respectively. Therefore, since the 
least squares routine did not provide acceptable start-up 
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values the identified parameters were set initally to zero 
for the first SCIDNT iteration. Further attempts to obtain 
phugoid start up values for SCIDNT from the least squares 
routine proved futile and its use was discontinued for the 
phugoid case. 

The version of SCIDNT utilized for this 
investigation would accept a maximum of 500 data samples per 
run. Since the phugoid period was approximately 100 sec a 
data burst recorded at 20 samples per second (SPS) would 
only include the first 25 sec of the maneuver. This time 
span was obviously acceptable for the short period mode but 
clearly unacceptable during a phugoid maneuver. An attempt 
was therefore, made to utilize data at 2 and 5 SPS to 
include as much of the phugoid maneuver as possible. This 
data was from the same maneuver with data tapes recorded at 
2 and 5 SPS. 

Although the 2 SPS data would provide a longer 
phugoid time span in theory, the program exhibited severe 
instabilities during the control input and therefore only 
the latter portion of the data burst was usable. It was 
strongly suspected that these instabilities were due to the 
unusually rough control inputs. 

The actual parameter identification procedure for 
the phugoid mode was somewhat different than that utilized 
during the short period investigation. The short period 
stability derivatives were obtained from a run whose initial 
conditions most closely approximated those of the phugoid 
run. These short period parameters were then fixed as 
constants during the phugoid iterations. Computer runs were 
made with the number of iterations and step cuts both 
initialized at ten. The program would then continue to 
iterate until the desired number of iterations/step cuts 
were completed. For both the 5 SPS and 2 SPS data 
convergence was obtained prior to executing the full 10 
iterations. Iterations conducted with phugoid data recorded 
at 2 SPS appeared to converge very quickly and began 
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executing step cuts until the allowable number was 
exceeded. As in the short period case the program executed 
one step cut more than requested. Figure 4. presents the 
change in the likelihood function during iteration using 
phugoid data recorded at 2 SPS. The program converged to its 
"best” solution by the third iteration at which time the 
step cuts were executed. 



Iteration 




Figure 4. Likelihood Function During Iteration 

Using 2 SPS 



Note that although the program had apparently 
converged at the third iteration the value of the likelihood 
function decreased until the fourth step cut. The function 
value then remained essentially constant throughout the 
remaining step cuts. 

The convergence characteristics using data recorded 
at 5 SPS were somewhat different. Figure 5. 
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Iteration 




Figure 5. Likelihood Function During Iteration 

Using 5 SPS 



It should be noted, however, that the scales are 

different for each plot. The program had essentially 

converged at the fifth iteration, at which time the program 

executed step cuts. Note that the likelihood function 

values decrease using 5 SPS data, and were essentially the 

same using the 2 SPS data. However, the time history matches 

for each set of data were significantly different. The angle 

of attack scale factor (K ) was not estimated since the 

' cc ' 

value of K from the short period data was utilized. After 
convergence was obtained for the parameter estimates a final 
set of iterations was performed to estimate values for 

Z , X , M , 0 and the three phugoid parameters 

o o c o 
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simultaneously. At this point, the phugoid time history fits 

were acceptable for the 5 SPS data and clearly unacceptable 
using the 2 SPS data. Further attempts to improve the 
estimates and time histories resulted in using a lot more 
computer time with no improvement. 

3 . Data Ana lysis 

Data from phugoid estimation runs using 2 SPS and 5 
SPS are presented in Tables VII and VIII 
respectfully. Phugoid time history fits are presented in 
Appendix D Figures 4 and 5 for 2 and 5 SPS respect fully. It 
should be noted that the scales of the phugoid time 
histories for 2 and 5 SPS are different. 



The data obtained using a sample rate of 2 SPS was 



totally 


unacceptable. For example, refering 


to 


Table 


VII, the 


value of H v/as -0.00C00017 while the 
u 


last 


step 


size was 


0.003. 
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Derivatives 



Z X X X M 

u OC u g u 



Estimate 


-.003 0 


-.758 


-47.9 


-1.7x10-7 0 


Standard 

deviation 


.00006 fixed 


.0542 


fixed 


.0002 fixed 


Parameter 
Step Size 


-.001 fixed 


0.444 


fixed 


0.003 fixed 


Eigenvalue 

Spread 




3.54 


X 108 




Likelihood 

Function 

Value 




■^8.4 


X 103 




Time 

History 

Matches 


Cr = 0.0048 
oc 


cr = 

u 


26.77 


(T = 0.00534 

g 




cr = 
e 


0.053 


cr = 2. 

n 
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Table VII. SCIDKT Estimates of Longitudinal Phugoid 
Stability Derivatives Using 2 SPS Data 



Although the eigenvalue spread and likelihood 
function change during iteration indicate no significant 
program difficulties the standard deviations of the state 
error vector are very large. As an approximation it can be 
said that the phugoid maneuver is essentially a varying 
airspeed and pitch attitude maneuver with a constant load 
factor. From Table VII it can be seen that the standard 
deviation for airspeed was 26.77 ft/sec~i and for pitch 
attitude was 0.053 rad, both of which were unacceptable. 

The time history plots constructed using the 2 SPS 
data showed no real agreement between the model and the 
airplane. In constructing these plots both the model and 
airplane were excited with the same control input and the 
airplane time history plots superimposed on the same 
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scale. It is obvious that there was simply no fit of any of 
the state vector. It should be noted that the plot of 
average horizontal stabilizer position although appearing to 
be very rough was a result of the plot scale only and was in 
essense relatively smooth. 

The parameter identification results using the 5 SPS 
data indicated it was possible to identify the phugoid mode 
parameters. Referring to Table VIII, it can be seen that the 
estimated derivatives using the 5 SPS data converged to 
appropriate (sign and magnitude) stability derivative 
values. 
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Estimate 


-.0002 
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-9.6x10-5 0 


Standard 

Deviation 


1.2x10-6 
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Table VIII. SCIDNT Estimates of Longitudinal Phugoid 
Stability Derivatives using 5 SPS Data 



The standard deviations of the parameter estimates were much 
smaller than the parameter, and the parameter step sizes 
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indicated that the estimates had converged to the "bast” 
estimates obtainable with the data. It should be noted that 
standard deviations for u and 6 were 5.53 ft sec~* and 0.014 
rad respectively, a significant reduction in standard 
deviation frcm those obtained with 2 SPS data 

The time history plots using 5 SPS data showed 
relatively close agreement between the airplane and model 
time histories. There was, however, some model instability 
noted during the first 30 seconds of the time histories. It 
was concluded that these instabilities were due to the 
relatively slow sampling rate during the initial control 
input. 

From SCIDNT the phugoid damping ratio ( Q ) and 

frequency (cO (rad sec~i)) were 0.0486 and 0.06 

n 

respectfully. Available F-14 data was obtained in the auto 

sweep mode with SAS on while SCIDNT data was obtained in 
auto sweep, SAS off. It was therefore not possible to make 
any meaningful comparisons. However, it was noted that 
between 10,000 ft and 35,000 ft the phugoid period varied 
between 106 sec and 72 sec. The SCIDNT frequency resulted in 
a period of 92.53 sec which when considered in light of 
available NPE data appeared reasonable. 

It was concluded that, within the scope of this 
research, the paramater identification program SCIDNT 
provided a reasonable estimate of the phugoid longitudinal 
parameters. 



C. SIKULTANEOUS SHORT PERIOD AND PHUGOID PARAMETER 
IDENTIFICATION 

1 . Airplane Test Condi tio ns 

Initial trimmed test conditions for the short period 
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and phugoid maneuvers are presented in Table IX. Prior to 
the control input the airplane was trimmed for steady state 
"hands off" flight [5]. Since ths short period input was 
applied after the phugoid input, the two test conditions are 
different. 



Altitude 


MACH 


CG 


Loading 


(Ft HPC) 




(% MAC) 


Phugoid Input 30 326 


.6217 


10.86 


B 


Short Period Input 30 703 


.5963 


10.86 


B 


Table IX. Short Period 


and Phugoid Test 


Conditions 



2* Method Of I nvestig a tion 

The airplane model equations used were those shown 

in Appendix B. The derivatives X , X^ , and X were fixed 

oc g 

as described in the previous section. 

The maneuver investigated in this section was 
designed specifically to provide the means of identifying 
the short period and phugoid parameters from one maneuver 
or, if you will, simultaneously. Initiation of the phugoid 
maneuver was accomplished in accordance with standard test 
techniques, [5], Upon establishing the airplane in the 
phugoid maneuver, a doublet input was applied to excite the 
short period node. The airplane was then allowed to fly 
through at least one phugoid cycle. 

Identification of the short period parameters was 
essentially routine and straightforward. The proceedures 
used in estimating the parameters were those described in 
the previous short period section. The computer runs were 
made with the number of iterations and step cuts both 
initialized at ten. Although this short period input was 
superimposed on the phugoid, convergence was obtained early 
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in program execution. Iterations conducted under these 
conditions converged very quickly and executed step cuts 
until the allowable number Was exceeded. Likelihood function 
change during the short period iterations was similar to 
that shown in Figure 2. 

Based on the results of the previous section it was 
concluded that data recorded at 5 SPS would most likely 
yield reasonable phugoid parameter estimates. The short 
period stability derivatives obtained from this run were 
fixed as constants during the phugoid iterations. Computer 
runs were then made with the number of iterations and step 
cuts both initialized at ten. 

Figure 6 presents the likelihood function value 
during iteration of the phugoid stability derivatives. 



Iteration 




Figure 6. Likelihood Function During Iteration 

For Phugoid 
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The likelihood function change for the phugoid maneuver 
displayed some characteristics heretofore not 

encountered. The characteristics shown are typical of the 
results obtained with the available data. At some time 
during the iteration sequence the likelihood function would 
exhibit a jump discontinuity, in this case it occurred 
between the first and second iterations. Program iteration 
then appeared to proceed in a normal manner until the tenth 
iteration at which time the program terminated. From Figure 
6 it was obvious that convergence had not been obtained at 
program termination. 

In an attempt to correct the iteration instabilities 

the phugoid stability derivatives Z , X , and M were 

u u u 

estimated simultaneously with Z , X , t5 , and 6 . Figure 7 

o o o o 

presents the likelihood function during a typical iteration 
for seven parameters. 
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I teration 




Figure 7. Likelihood Function During Iteration For 

Phugoid (Stability Derivatives, Z , X , I'l , and 6 ) 

o o o o 

It can be seen that the likelihood function value converged 
rapidly to a steady state value with no apparent 
discontinuities or nonlinearaties. 

3 . Data A nal ysis 

Data from a typical short period run superimposed on 
the phugoid are presented in Table X. A typical time history 
fit is presented in Appendix D, Figure 6. 
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Derivatives 
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Estimate 


-.4285 


1. 048 


-1 . 54 


-.665 


0.19 4.464 


Standard 

Deviation 


0.0049 


0.0074 


0.014 


0.011 


0.0082 0.0312 


Parameter 
Step Size 


0.0149 


0.0119 


-.0166 


0.0355 


0.0173 -.1815 


Eigenvalue 

Spread 






1.25 


X 105 




Likelihood 

Function 

Value 






-4.36 X 103 




Time 

History 

Matches 


cr = 0. 
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00188 


O' = 
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0.8822 


(T = 0.0036 

g 






0~ = 
0 


0.0048 


(T = 1. 
n 
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Table X. SCIDKT Estimates of Longitudinal Short Period 

Stability Derivatives 



The short period parameter estimates displayed 
characteristics similar to those observed when estimating 
the short period parameters from a short period input 
only. The time history plots showed good agreement between 
the model and actual airplane data. It was noted that 
airspeed was not constant since the short period was 
superimposed on the phugoid. This varying airspeed did not 
however, have any noticable effect on the determination of 
the short period parameters. 

The results for the phugoid portion of the maneuver 
were at times in disagreement with theory. From theory it 
can be shown that the speed damping derivative, which acts 
along the negative X axis increases with speed and is 
therefore negative [6]. During the iteration to determine 
the phugoid stability derivatives it was noted that the 
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speed damping derivative (X ) , became positive on the 

u 

seventh iteration and remained positive until program 
termination. Table XI. 



Derivatives Z X 

u oc 



Estimate -.0002 0 



Standard 1.1x10~* fixed 
Deviation 



Parameter 5.38x10~^ fixed 
Step Size 



Eigenval ue 
Spread 



Likelihood 

Function 

Value 



X X M X^ 

u g u *^e 

.005 -43.85 -7.8x10-6 0 

.0002 fixed 3.5x10-6 fixed 

4.8x10-5 fixed 1.2x10-5 fixed 

6.2 X 1012 

-1.005 X 10-» 



Time 

History 

Matches 



(T = 0.004 cr = 1.97 Cr = 0.0047 

oc u g 



O' = 0.031 cr = 1.48 
0 n 



z 

Table XI. SCIDNT Estimates of Phugoid Stability 

Derivatives (Positive X ) 

u 



The time history plots were in reasonable 

agreement. Appendix D Figure 7. The positive value of X was 

u 

not however considered reasonable since this value would 

imply decreasing drag with increasing airspeed. The X term 

u 

was however negative through the sixth program iteration. As 

mentioned previously the likelihood function during 
iteration showed definite signs of instability, it was 
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therefore suspected that the speed damping sign change was 
due to this program instability. The program was 
therefore, run again and terminated at the sixth iteration 
thereby producing a negative speed damping term. Table XII. 



Derivatives Z 
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oc 
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u 




Estimate -.00008 


0 


—.00032 


-43.85 


0.00017 


0 


Standard 2.6x10~* 

Deviation 


fixed 


0.00047 


fixed 


6.7x10-6 


fixed 


Parameter 4.0x10~s 
Step Size 


fixed 


-.0059 


fixed 


0.000058 
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Eigenvalue 

Spread 
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X 10 »i 






Likelihood 

Function 

Value 




-7.68 


1 X 103 
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Table XII. SCIDNT Estimates of Phugoid Stability 

Derivatives (Negative X ) 

u 



With the iterations thus restricted to six the value for X 

u 

was maintained at a negative value. Time history plots for 

the negative X case are presented in Appendix D, Figure 
u 

8. It was noted that there was no appreciable improvement in 

the time history plots. In fact from Tables XI and XII it 
can be seen that the standard deviation of the state error 
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case. This 



vectors was degraded for the negative X 

u 

condition in fact is however, not very startling when it is 

realized that the program had demonstrated instabilities and 
had not converged. 

In an effort to reconcile these inconsistencies and 

lack of convergence the three phugoid parameters were 

estimated along with Z , X , H , and 6 . This technique 

o o o o 

yielded reasonable stability derivatives and acceptable time 

history plots. Stability derivative data provided by the 
seven parameter identification are shown in Table XIII, time 
history plots are presented in Appendix D, Figure 9. 
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Estimate - 


.00019 


0 


-.0029 -43.85 


-.000141 


0 


Standard 1 

Deviation 


.hx10-6 


fixed 


.00026 fixed 


5.2x10-5 


fixed 


Parameter 1 
Step Size 


.23x10-7 


fixed 


6.3x10-5 fixed 


-4.02x10-6 
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Eigenvalue 

Spread 






3.7 X 1013 






Likelihood 

Function 

Value 






-4.23 X 103 






Time 

History 

Hatches 


• 

o 

li 


0045 


C 

It 

• 


<r =0.00523 

q 






or 

e 


= 0.031 O' = 1. 

n 


86 





z 



Table XIII. SCIDNT Estimates of Phugoid Stability 

Derivatives (Parameters Estimated With Z , X , M , and 6 ) 

o o o o 
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The parameter standard deviations were small and 
within engineering tolerance. The parameter step sizes were 
quite small and showed further evidence of solution 
convergence. The standard deviation of the calculated model 
response with respect to the measured airplane response was 
acceptable. The time history plots provided further 
indication of parameter estimate validity. 

The problems associated with the phugoid using the 
simultaneous technique could have been a result of program 
deficiencies or airplane anomalies or a comnination of 
both. The flight conditions at which the data was obtained 
were areas in which the phugoid mode was lightly damped or 
neutrally damped. In addition the airplane wing was set to 
the auto sweep node. It was felt that either of these 
conditions could have caused the program instabilities 
noted. Attempts to obtain new data were not possible due to 
airplane non-availability as mentioned earlier. 

It was concluded that extraction of the short period 
stability parameters in the presence of the phugoid was 
easy, straight forward, and yielded results similar to those 
obtained from pure short period data. Estimation of the 
phugoid stability parameters was possible when they were 
estimated in conjunction with Z , X , K , and 6 . 



It 

obtained at 
moderately 
the present 



o o o o 

was recommended that a new set of data be 
a flight condition where the phugoid is at least 
damped and that this data be analyzed to resolve 
anomalies. 
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III. CONCLUSIONS 



A. GENEBAL 

The maximum likelihood identification program SCIDNT 
showed great promise in obtaining the short period and 
phugoid stability derivatives from one maneuver, 
"simultaneously". 

B. SPECIFIC 

Within the scope of this research, the parameter 
identification program SCIDNT was easy and straight-f crward 
to use in determining the short period mode. 

Meaningful parameter estimates were obtainable only when 
started from an initial steady state condition. 

Within the scope of this research, the parameter 
identification program SCIDNT provided a reasonable estimate 
of the phugoid longitudinal parameters. 

Extraction of the short period stability parameters in 
the presence of the phugoid was easy, straight-forward, and 
yielded results similar to those obtained from pure short 
period data. 

Estimation of the phugoid stability parameters was 
possible when they were estimated in conjunction with 

Z , X , M , and 0 . 
o o o o 
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IV. 



RECOMMENDATIONS 



It is recommended that a new set of data be obtained at 
a flight condition where the phugoid is at least moderately 
damped and that these data be analyzed to resolve the 
present anomalies. 
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APPENDIX A 



TEST AIRCBAFT DESCRIPTION 



The F-14 airplane is a supersonic, two-place, 
twin-engine, swing-wing, air-superiority fighter designed 
and manufactured by Grumman Aerospace Corporation. In 
addition to its primary fighter role, carrying missiles and 
an internal 20-millimeter gun, the airplane is designed for 
fleet air defense and ground attack missions. Photographs of 
the test airplane are shown in Figures 1 and 2. A three-view 
drawing is presented in Figure 3. 

The forward fuselage, containing the crew and electronic 
equipment, projects forward from the mid fuselage and wing 
glove. Outboard pivots in the highly swept wing glove 
support the movable wing panels, which incorporate integral 
fuel cells and full-span leading edge slats and trailing 
edge flaps for supplemental lift control. In flight, the 
wings may be varied in sweep, area, camber, and aspect ratio 
by selection of any leading-edge sweep angle between 20° and 
680. Wing sweep can be automatically or manually controlled 
to optimize performance and thereby enhance aircraft 
versatility. Separate variable-geometry air inlets, offset 
from the fuselage in the glove, direct primary airflow to 
two TF30-P-412A, dual-axial-compressor, turbo-fan engines 
equipped with afterburners for thrust augmentation. The 
displaced engine nacelles extend rearward to the tail 
section, supporting the twin vertical tails, horizontal 
tails, and ventral fins. The middle and aft fuselage, which 
contains the main fuel cells, tapers off in depth to the 
rear where it accommodates the speed brake surface and 
arresting hook. Retractable vanes in the glove leading edges 
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extend to supplement lift and compensate for changes in the 
aircraft aerodynamic center. All control surfaces are 
positioned by irreversible hydraulic actuators to provide 
desired control effectiveness throughout the flight 
envelope. Stability augmentation features in the flight 
control system enhance flight characteristics and provide a 
more stable and maneuverable weapons delivery platform. The 
tricycle-'ty pe, forward- retracting landing gear is designed 
for nose gear catapult launch and carrier landings. Missiles 
and external stores are carried from eight hardpoint 
stations on the center fuselage between the nacelles and 
under the nacelles and wing glove; no stores are carried on 
the movable portion of the wing. The fuel system 
incorporates both inflight and single-point ground refueling 
capabilities. 

A more complete description of the F-14A airplane is 
given in Ref. 8. 
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SIDE VIEW 



F-14A AIRPLANE 



Figure 2 
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APPENDIX B 



DERIVATION OF THE AIRPLANE MODEL AND MEASUREMENT EQUATIONS 

A detailed derivation of the complete aircraft equations 
of motion is shown in Ref. 6. A summary of the assumptions 
made during that derivation is shown below for 
completeness. These assumptions lead to a set of force and 
moment equations. 



* Summary of Derivation Assumptions * 

1. The airframe is a rigid body 

2. The earth is fixed in space 

3. The earth's atmosphere is fixed with respect to the 
earth 

4. The mass of the airplane remained constant for the 
duration of any particular dynamic analysis 

5. The xz plane is a plane of symmetry 

6. The disturbances from the steady flight condition is 
small enough so that the products and squares of the changes 
in velocities are negligible in comparison with the changes 
themselves 

7. The disturbance angles are small enough such that the 
sines of these angles may be set equal to the angles and the 
cosines set equal to one. Products of these angles are also 
approximately zero and can be neglected 

8. Since the disturbances were small, the change in air 
density encountered by the airplane during any disturbance 
can be considered zero 
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9. During the steady flight condition, the airplane is 

flying with wings level and with all components of velocity 

zero except U and W 
o o 

10. Body axis system is assumed 



Dividing the force equations by the mass ra and the 
moment equations by the appropriate moments of inertia 
yields terms of the form: 

J_4X_u_ - and 1 JL r. ( B1 ) 

m "p 3^^ 

XX 

Replacing the first term by X and the second term by L 

u r 

simplifies the notation. These quantities are called 

"dimensional stability derivatives". Using this short hand 

notation the equations of motion applicable to the 

longitudinal problem can be written as follows: (^, complete 

set of equations is presented in Ref. 5) 



• 

u 


+ w 


q 


+ ge 


cos 0 


= X 


u 


+ 




• 

u 


+ X q + 


x.g 






o 




o 




u 




u 




g 


g 


+ 


X w 


+ 


X.w + 


CP 


+ 




b 




+ 


Xtf o 






w 




V 


e e 




e 




e 




e e 






x<y 


6 


+ X^ 


b + 




6 


+ 


X 


• 


• 

a 






F 


F 


F F 


°B 


B 






B B 






+ 


^'6 


• • 


cos c 


T u 


+ 


cos J 


X 


T(5 6 








®B 


B 


o 


u 








> 


^RPM 


RPM 



80 
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In order to simplify these equations it was assumed that 
all disturbances were small such that the higher order terms 
were negligible. Also it was assumed that the flaps and 
speedbrakes were retracted and that thrust effects were 
negligible. With these further assumptions the w equation 
became: 

w - 0 g + g6 sin0 = Zu+Zg + Zw + Z^( 5 > (B5) 

o o u q w ^e e 

It was further assumed that the flow was 

guasisteady. Because of this all derivatives with respect to 

the rates of changes of velocities were omitted with the 

exception of those with respect to w which were retained to 

account for the effect on the horizontal tail of the 

downwash from the wing. If a wing moving through an air mass 

with an angle of attack aC and a steady forward velocity U 

o o 
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is suddenly given a vertical velocity change w without 
changing the forward velocity, the angle of attack 
perturbation is as shown in Figure 1. 





Figure 1. Angle of Attack Due to Vertical Velocity 

Perturbation 



Thus changes in angle of attack, , are proportional to 
the vertical velocity perturbation w. Thus we can write. 



w + W 

o 

tan(OC + oc ) , 

o *VU2+ (w + W ) 2 



and by the snail angle assumption 
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oc = 
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and oc 
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Taking the derivative yields. 



w = U OC 
o 



also it can be shown that. 
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and substituting equation B 8 yields. 



3Z 1 aZ 

= Limit — 

0w oc — ►o U ^ 
o 






also if we redefine the Z dimensional derivatives 
typically. 
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ccucc uuu giig 
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we can now substitute into equation B 3 which 
equation , 



^=ZOC + Zu + Zq e sin 0 + Z^ ^5 

cC u q U o e e 

o 



where Z = 1 + Z 

q g 



A similar expression was obtained for u 
simplifying assumptions and the expressions. 



SX 1 aX 1 

= Limit — - = X 

0w oC— +>o U OC u C5C 

o o 



and the expression for u was therefore written. 



u = X CC+Xu+X^q -gG cosG + X ^5 

OC u q o e e 



( BIO ) 
such that, 

( B11 ) 
yields the 

( B12 ) 

using the 
( 313 ) 

( B14 ) 
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where x' 

g 



= (X - cc U ) 
g o o 

» 

Finally the expression for g followed directly as, 

g = M OC +M u + a q +M. <S ( B15 ) 

^ u g *^6 e 

The model equations for the longitudinal case, including 
an angle of attack due to gusts term, were written in matrix 
form as. 



X = Px + Gu + Pw 



( B16 ) 
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The measurement equations written in matrix form were. 



y = Hx + Du + V 



( B17 ) 
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Note: Selected stability derivatives are primed since 
they differ from the notation used in Ref. 6, however they 
are exactly equivalent to the derivatives used in SCIDNT. 
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APPENDIX C 



AIRPLANE AND GROUND DATA SYSTEMS 



A. AIRPLANE DATA SYSTEM 
1 , G eneral 

The demonstration F-14A ship number 8 BU. NO. 157937 
was the test airplaive utilized. A detailed description of 
the airplane instrumentation system is contained in Refs. 7 
and 9. The airplane was modified to be compatible with the 
Naval Air Test Center Real-Time Telemetry Processing System 
(RIPS) . The instrumentation systems and measurements 
selected for use in the test airplane were primarily 
designed to satisfy the Flight Test requirements for the 
F-14A Performance Evaluation and Demonstration. The test 
airplane was not entirely a production configured. F-14A in 
that some miscellaneous cockpit components not required for 
flight testing were removed. The production nose radome was 
replaced by a fiberglass and metal radome fitted for an 
instrumentation noseboom assembly. Additionally, some 

production equipment was reworked in order to accomodate 
various instrumentation subassemblies. 

Nose Boom Inst rume nt ation 

The nose boom assembly included Angle of Attack and 
Angle of Side Slip Vanes and a Pitot / Static Head. 
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3 • Co ckpit Inst rumentat ion 

The Forward Cockpit instrumentation was primarily 
pilot oriented and included; Angle of Attack and Sideslip 
Indicators, Airspeed, Altimeter, Machmeter, Center of 
Gravity Normal Acceleration and an Instrumented Control 
Stick. The Aft Cockpit Instrumentation included left and 
right engine EPM (N1), Main Fuel Flow, A/B Fuel 
Flow, Outside Air Temperature and the Instrumentation Master 
Control Console. 

4 . Compa rtment Instru me ntat i on 

The instrumentation subcomponents were mounted 
throughout the airplane in various available locations 
including, the Nose, Forward Compartments, and Equipment 
Compartments. The bulk of the instrumentation equipment was 
located in the forward equipment compartments. 

5 . Sxs t € m Descri p ti on 

The instrumentation system consisted of the 
following subsystems: 

o PCM Signal Conditioning 
o PCM 
o Timing 

o Magnetic Tape Recording 
o Telemetry 

o Power and Distribution 
o Ground Support Equipment 

The system had the capacity to record on magnetic 
tape and to telemeter to the ground-data-reduction-ccmplex 
all the PCM data for real-time analysis. The PCM Signal 
Conditioning Subsystem consisted of active pre-sample filter 
chassis and seme special conditioning required for certain 
measurements. The pre-sample filters had cut-off frequencies 
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of 5^ 10, 20, 40, and 80 Hz, The PCa subsystem consisted of 
an analog multiplexer, analog to digital converter, digital 
multiplexer, and the required control logic. The subsystem 
had the capacity of accepting up to 332 analog parameters at 
frequency response from dc to 40 Hz. Discrete data (On/Off 
signals) and a limited amount of digital data and time words 
were time-division multiplexed into the output data 
stream. The output data was presented as a succession of 
twelve-bit (including parity) parallel data words in NRZM 
format which were recorded on the digital tape. A serial 
output in NRZL format for entry into the ”L" Band telemetry 
subsystem was provided. The PCM format is contained in Ref. 
10. An IRIG-c Time Code Generator was installed to provide a 
common time correlation for magnetic tape and telemetry 
data. The Magnetic Tape Recording Subsystem consisted of a 
PCM Recorder which had the capability to record the PCM data 
and Time Code information. The Digital Recorder contained a 
sixteen tracks IRIG compatible head capable of storing one 
hour of digital data on one inch tape at a tape speed of 30 
IPS. The data was recorded in parallel using saturation 
recording techniques. This portion was capable of bit 
densities of 1000 bits per inch. The i'elemetry Subsystem 
consisted of a signal combining unit, an '•!” Band 
Transmitter and two blade antenna. The operating frequencies 
could be tuned to coinside with those in use at various test 
facilities. The entire Flight Test Instrumentation System 
was powered from the aircraft primary (generating) 
system. Should either the left or right engine or generator 
fail, the instrumentation system would continue until such 
time as both generators failed. The distribution system 
consisted of a 28V dc transformer/rectifier (T/R) , and a 
power/control distribution box, which provided electrical 
system protection and isolation from the aircraft power 
systems through individual circuit breakers. The equipment 
required to perform checks on the installed instrumentation 
systems was configured in mobile test carts. These carts 
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required only 60 cycle power for operation and provided the 
operator with a tool for conducting pre-flight and post 
flight tests on the tape- recording system. 

B. GROUND DATA SYSTEM 

The Naval Air Test Center, Patuxent River, Maryland, 
employs a Real-Time Telemetry Processing System (RIPS) as 
universal support for testing on several types of aircraft 
and weapons systems, Ref. 10. The computer system, designed 
by Xerox, consists of a Sigma 9 general purpose computer as 
the central processor supported by three Sigma 3 computers 
performing the data input and output processing. RTPS 
provides the capability to process two complete and 
independent tests (data "streams”) simultaneously from 
either tape or telemetry sources. These "streams" can 
support totally unrelated and asynchronous test flights or 
can support setup for one flight while processing another 
flight. 

Each Data Channel accepts up to 512 channels of PCM data 
and 32 parameters of FM data. This data is processed 
(individually or mixed) by a Sigma 3 and preprocessing 
Programmed Algorithm Unit (PAU) . Up to 50,000 values of 
Engineering Units (EU) data per second can be fed to the 
Sigma 9 while simultaneously being stored on digital tapes 
for later playback. Two large screen interactive graphic 
display stations, that provide real-time "hands on" data 
analysis and presentation, are controlled by the Sigma 3 
Display Processor. Batch jobs can be processed when the 
real-time and non real-time operational modes are not 
utilizing the resources of the system to 100 percent. The 
batch job is defered if additional demands are made by 
either non real-time or real-time tasks during batch job 
processing. 

The specific flight test data utilized by the parameter 
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t 



identification 
pitch attitude. 



program were, V , ALPHA PEL, pitch rate, 

T 

NZ NET, and average horizontal stabilizer 



position. The raw data (digital) was used to generate an EU 
tape for input to the paraaeter identification program, in 
generating the EU tape corrections were made to the measured 
values as required in accordance with Ref 11. For example 
corrections tc angle of attack included nose boom bending 
moments, fuselage bending moments and pitch. 
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APPENDIX D 



TIME HISTORY FITS 





Short Period Time History Using Doublet Input 



Figure 1 
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